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EXECUTIVE SUMMARY • 

The severe drought experienced In Central Florida throughout most of 
1981 brought Increased attention to the long recognized problem of fresh 
water availability across the State. Much of the state's fresh water supply 
Is obtained from groundwater storage and natural lakes. Many of these lakes 
reached record and near-record low levels during late sunnier of 1981 and 
forced the Imposition of the restricted water use practice for Individual 
households and municipal supplies as well as agricultural irrigation. 

Management of water resources In Florida Is the responsibility of 
Florida's five Water Management Districts. Especially during drought condi- 
tions. a close watch Is maintained over all water supplies. This Is most 
comtionly done through monitoring water stage, the water height above sea 
level. Although water stage Is an Important Indicator, water storage volume 
Is a more critical parameter. Volume can only be obtained If accurate data 
exists as to the lake contours in conjunction with water stage. Such data 
rarely exists because of the extensive time and cost associated with conven- 
tional transect survey methods of collecting such data and the fact that they 
are subject to change over extended periods of time. 

The University of Florida. In cooperation with the St. Johns River 
Water Management District, and Kennedy Space Center developed a technique 
using Landsat data for estimating available water storage volume and applied 
It successfully to Lake Washington and Lake Harris in central Florida. A 
number of Landsat scenes including the lakes of Interest were selected to cor- 
respond with a wide range of lake stages as measured over the past nine years. 
Lake surface area was then measured from the Landsat data, and when properly 
averaged, was used with the change In lake level to estimate the change in 
lake volume. Thus, a change in lake stage can be directly correlated with 
available water volume. 


Eight cloud dfltos woro chosen for Lake Washington, Florida. Tho 
lake stages on those dates were recorded. The dates were rearranged in 
order of Increasing lake stage. The water surface was measured from Uandsat 

along with the lake stage that ranged from 10.60 to 15.85 ft. (msl). The re- 

suits Indicated that the water surface varied from 2,537 acres at the stage 10.60 
ft. to 2,850 acres at the stage 15.85 ft. The ground truth measurement of water 
surface around, the stage 15.85 ft was 2848 acres which was very close to 2850 
acres, the area estimated by the Landsat system at the stage 15.85 -ft. The 
difference was negligible. This Implies that the technique developed In this 
study Is quite applicable to measure the water surface area of the lake. The 
lake volume between two stages was computed from the average area between the 
stages multiplied by the stage increment. The lake volume between the stages 
10.60 and 15.85 ft Is about 14,352 acre-ft. The lake volume ranged from 26 to 
28 acre-ft. for each 0.01 ft. Increment in lake stage. 

Four cloud free dates were chosen for Lake Harris, Florida. The similar 
techniques as used In Lake Washington were applied to Lake Harris for measuring 
the lake water surface area. The results Indicated that the water surface 
varied from 17,430 acres at the stage 62.38 ft (msl) to 18,657 acres at the stage 
63.30 ft. The lake volume between the stages 62.38 ft. and 63.30 ft. was about 
16,682 acre-ft. The lake volume ranged from 176 to 183 acre-ft. for each 
0.01 ft. change In lake stage. 

The rates of change In lake volumes In relation to lake stage Increments 
are quite stabilized In both lakes. This Implies that the lake volume at 
other stages could be roughly estimated based on the lake volume Increment re- 
lation obtained In this study. 

The key element In this approach, I.e. the lake surface measurement from 
Landsat was accomplished by density slicing In the near infrared-band 7. For 
lakes having shallow depths and marshy shares as represented here, such mea- 
”'SurementS"’ai'e"h1^>'ly'- sensitive to the upper limits chosen for band . The 




upper limit of 35 was selected In this study. 

The technique developed 1'n this study can be applied two ways, First* 
where the historical stage records are available, the historical Landsat data 
can be us'sd to establish the relationship between lake volume and lake stage. 
In the second case, where the historical stage records are not available, the 
historical Landsat data can be used to estimate the historical lake stage 
after the lake volume and lake stage Information becomes available In the 
future. 

General Information on Landsat remote sensing system Is presented In 
Appendix A In this report. Computerized methods used for bounded area and 
weighted sub-area computation are presented In Appendix B. 


m4» 

introduction 

Florida Uglslaturo passed the Florida Water Resources Act of .1972,. which 
established five Water Managomont Districts in the. state of Florida and re- 
quired that each district should develop a Water Use and Supply Development 
Plan that would take Into account all factors Involved up through the year 
2020. To moot the mandate of this legislation the districts have, to use every 
available technology and are looking toward National Aeronautics .and Space 
Administration (NASA) remote sensing and Automatic Data Process (ADP). techniques 
for major help In solving the complex problems. Even though NASA Is presently 
Involved In considerable research In the area of water resources, much of this 
work cannot be directly applied to the problems of the area. For Instance, 
the geology of the district with Its sandy soil, shallow lakes, and large flat 
marshland areas make the storage of adequate surface water difficult. These 
water conditions not only can produce severe water shortages during the winter 
dry season which coincides with a period of high agricultural water demand, 
but also can cause an acute problem of flooding during wet season. For Instance 
the severe drought experienced In Central Florida throughout most of 1981 and 
the beginning of 1982 brought increased attention to the long recognized. problem 
of fresh water availability across the State. Much of the state's fresh water 
supply is obtained from groundwater storage and natural lakes. Many of these 
lakes reached record and near-record low levels during late summer of 1981- and 
forced the imposition of the restricted water use practice for individual 
households and municipal supplies as well as agricultural irrigation. 

A typical basin in central Florida Is the St. Johns River basin (as shown 
in Figure 1) which encompasses approximately 11,430 square miles. The river 
rises in poorly defined marshland west of Fort Pierce and Vero Beach in east 
central Florida. It flows northward, roughly parallel to Florida's east coast 
for 285 miles before merging into the Atlantic Ocean just north of Jacksonville 
Management of water resources in St. Johns River basin is the responsibility 
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of St. Johns River Water Management Dls-trlct (SJRWMD). A close watch Is 
maintained over all water supplies, especially during drought conditions. -This 
Is most commonly done through monitoring water stage, the water height above 
sea. level. Although water stage Is an Important Indicator, however, water 
storage volume Is a more critical parameter for proclso assessment of available 
water. Volume can only be obtained If accurate data exists os to the lake 
contours In conjunction with water stage. Such data rarely exists because of the 
extensive time and cost associated with conventional transect survey methods of 
collecting such data and the fact that they are subject to change over extended 
periods of time. Therefore, the management of water within the basin poses 
many unique and difficult problems that warrant the use of recently developed 
Landsat remote sensing techniques to help tK’. district make realistic decisions 
concerning water resources planning .and management. 

The University of Florida, In cooperation with the St. Johns River Water 
Management District, and NASA, Kennedy Space Center Initiated a joint project 
to deal, with the district-wide Water Resources Investigations and Management 
using Landsat Data. The Phase I of this study Includes the lake volume determi- 
nation. Recently^ applications of Landsat data to determine the lake volume In 
Lake Okeechobee have been done by Shih (1980), and Gervin and Shih (1981). The 
basic technique used In their study was to analyze the digital data from the 
Landsat earth-orbiting satellite by using the General Electric's multlspectral 
Image analyzer, the Image 100. In typical classification analysis, an area 
where land cover Is known from ground observations Is located as a training site 
on the Image and its spectral characteristics are measured. Then all areas which 
have similar spectral characteristics are Identified and assigned a color and a 
number code. The similar spectral characteristics are ,'eferred to as the sig- 
nature for that land cover type and the location of the color-coded pixels 
possessing that signature Is displayed on the CRT screen as a theme. Once the 
signature Is finalized, the pixels (picture elements) In the satellite scene 


which posses that signature (theme) can be counted. Thus, the area of that 
signature Is obaalned. A more, detailed description of the Landsat data 
analysis was presented by Gervin and Shih (1981). The general Information 
on Landsat remote sensing systems 1s presented as Appendix A In this re- 
port. However, the previous studies reqiil red some ground truths of the vegetation 
classification before a proper area could be estimated from each signature. 

In other words, the sole band which can be used to identify the land-water 
boundaries directly without the ground truth Information was not emphasized 
well previously. Therefore, In this study, the Image enhancement of sole band 
with density slicing technique was used to determine the lake water surface 
area. 

In developing a new technique to be used for Landsat remote sensing appH- 
cation, a groundtruth calibration is an Important procedure to demonstratp '•‘le 
applicability of the new technique. Thus, computerized methods - t .u 
area and weighted sub-area computation were also presented in this study. 

The general objective of this study was to demonstrate the application of 
Landsat remote sensing techniques for determining the lake volume of Lake 
Washington (Figure 2), and Lake Harris (Figure 3). The specific objectives were: 

1) to introduce the general information on Landsat remote sensing systems; 

2) to introduce computerized methods used for surface area and sub-area 
computation; 

3) to develop a new method used to compute water surface area from the 
Landsat data; and 

4) to discuss the applicability of the newly developed method for lake volume 
determination using Landsat data. 


MATERIALS AND METHODS 

Two aspects are Involved In this section; remote sensing techn1<^uos 
and groundtruth calibration. . 

Remote Sensing Techniques 

Landsat Remote Sensing Systems ; The general Information on Landsat 
remote sensing system Is presented as an Appendix A In this report. Nine 
sections comprise Appendix A: 1) introduction to Landsat Remote Sensing 

System; 2) Ssnsor, Band Designation, Wavelength, and Resolution which Include 
the Return-Beam Vidicon (RBV) Camera, Multlspectral Scanner System (MSS), and 
Thematic Mapper (TM); 3) Ground Receiving Stations; 4) General Electric 
“Image lOD" System; 5) United States Distribution Center; 6) The Northeast 
Regional Data Center, Florida; 7) University of Florida, IFAS Remote Sensing 
Facilities; 8) Landsat Newsletters; and 9) Landsat Update Information. 

As mentioned In Appendix A, Landsat passed over the same swath (185 km 
wide) of the ground every 18 days during the period from January 23, 1972 to 
January 22, 1975, after that the coverage period was 9 days except during the 
period between January 6 ana March 5, 1978 when the operation of Landsat 1 
was ended and the Landsat 3 was not launched as yet. Among the sensors of 
Landsats 1, 2, and 3, Multlspectral Scanner Systems (MSS) have provided the re- 
mote sensing community with a preview of what- to expect from an operational 
remote sensing satellite. Certainly, the MSS system Is going to be replaced 
by the Thematic Mapper system on Landsat 4 system. The MSS systems on Landsats 

1, 2, and 3 are operated on bands 4, 5, 6, and 7 and another band 8 Is also operated 
In Landsat 3. The Thematic Mapper systems on Landsat 4 are operated on bands 1, 

2, 3, 4, 5, 6, and 7. The Landsat 4 Is planned to be launched during the summer 

of 1982. 

Before 1979, the Landsat data was stored on computer compatible tape (CCT) 
which was produced at the Goddard Space Flight Center Image Processing Facility 
from wideband videotape recordings received by stations of the Satellite Tracking 


and Data Network. The Image data are reformated on h.1gh density digital 
tape. Since 1979, a new system has been operated to shorten the processing 
time. The user can have the tapes from the User Services Section. U.S. 

Geological Survey, Earth Resource Observation System Data Center* Sioux 
Falls. South Dakota, 57198 and the organizations mentioned In the 
Appendix A. 

As the Landsats have been In operation for more than nine years, a large 
amount of data has been accumulated on the Lake Washington and Lake Harris. 

In order to save Image analyzing time, selection of useful tapes Is necessary. 

The procedure to pick up useful tapes can be summarized as follows: 

1. Calculate the date of which Landsats flew over the Lake Washington. 

and Lake Harris. 

2. Pick up the cloud free days from the. dates obtained In. the. previous 
step by consulting meteorology. 

3. Pick up the various dates, when the water level Is distributed from 
low to high levels. In this step, for a given observation date the 
lake water stage Is required. 

Using Landsat Data for Water Surface Area Measurements. As mentioned in the 
Appendix A, the characteristics of band 7 provide best penetration through 
haze and light clouds and it also emphasizes live vegetation and land-water 
boundaries. For these reasons, band 7 was chosen in this study for the lake 
water surface area measurement. 

The Landsat high density digital tapes of the chosen dates for the Lake 
Washington and Lake Harris were analyzed on General Electric's multlspectral 
Image analyzer, the Image 100. The Image 100 consists of two tape drives, for 
Inputting and outputting data; a color CRT screen which displays (512)^ pixels 
of digital data; a memory for storing and refreshing the displayed Image, and a 
battery of programs capable of measuring, manipulating and highlighting the 
satellite data. Because the spectral characteristics of a water body, vegetations 
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and ground are different, each characteristic can be Identlfleo and assigned, 
a color and a number code. The spectral characteristics are referred to as the 
signature for the land cover type and the location of those color-coded pixels 
possessing that signature Is displayed on the CRT screen (or mapped on paper) 
as a theme. Once the signature Is finalized, the pixels In the satellite scene 
which possess that signature (theme) can be counted. In this way, maps and 
tabulations of areas that have similar spectral properties can be obtained. 

The key element 1n this approach, lake surface measurement from Landsat was 
accomplished by density slicing In the near Infrared band 7. The density slicing 
Is a procedure used to carry out the enhancement of the remote sening data. 

The density slicing Is a digital approach which has the advantage of flexibility. 
Any level or levels of grey may be selected on a photographic or electrical 
Imprint, with output via a line printer, on a flat-bed or drum plotter, or directly 
onto a photographic film writer, CRT screen. Quantization noise or spurious con- 
touring may result If the slices are not chosen carefully. For lakes having 
shallow depths and marshy shores as represented here, such measurements are 
highly sensitive to the upper limits chosen for band 7. The upper limit of 25 
was selected In this study. 

Lake Volume Computation ; After lake water surface area was estimated by 
Landsat data, the following equation was used to calculate the Increase In lake 
volume associated with a small Increase In lake stage: 

Vl “ Vg + (1) 

where » lake storage volume at lake stage s + 1 In acre ft; 

Vj ■ lake storage volume at lake stage s. In acre-ft; 

^s+1 “ surface area at lake stage s + 1 in acres as 

measured by Landsat data; 

Ag » lake surface area at lake stage s as measured by 
Landsat data, in acres; and 

Ah » the range between lake stages. 
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Because of the Initial lake volume at s«l Is unknown, the Vj»C Is used, 

where C is a constant. 

Ground Truth Measurement . .. 

A contour map provided by the St. .Johns. River Water Management District 
was used to compute the surface area. A computerized method as presented 1n 
Appendix B was used to compute the ground truth area measurement. In order to 
explore further potential applications of computerized methods for bounded 
area and weighted sub-area computations the detailed computer programs for each 
method are also presented as an Appendix B_1n this report. 
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RESULTS AND DISCUSSION 

Dates Selection 

As presented In the previous section* three step procedures were used to 
select cloud free dates associated with a wide range of lake stage fluctuations 
for areas of Lake Washington and Lake Harris, Florida. The results for selected 
eight cloud free dates of 9/6/72, 11/29/73, 2/27/74, 3/17/74, 6/15/74, 10/19/74, 
2/9/76 and 4/11/76 were used for Lake Washington and those used for Lake Harris 
were taken from four cloud free dates of 9/6/72, 8/31/75, 2/14/75 and 1/22/76. 

The lake stages on those dates were also recorded. 

Area Measurement 

Lake Washington ; The dates were rearranged as shown in Table 1 In the In- 
creasing order of the lake stage. The Lake Washington water surface area at 
those eight dates was estimated from the Image 100 after density slicing of the 
band 7 was done and a theme for the lake water surface was assigned. A typical 
color-coded classification result of water surface area for the Lake Washington 
on date 10/19/74 is shown In Figure 4. The results of water surface area esti- 
mated by Landsat data on those eight dates are shown in Table 1. The surface area 
varied from 2,537 acres at the stage of 10.60 ft above mean sea level to 2,850 
acres at the stage 15.85 ft. The stages are plotted against the water surface 
area as shown on Figure 5 which can be used to estimate the surface area associated 
with any stage between the measurements. For instance, the surface area Is 2816 
acres when the stage is 14.6 ft. 

The ground truth measurement was made from the contour map provided by 
the District using the computerized method as given in the attached Appendix B. 

The computed value was 2,848 acres around the stage 15.85 ft. This is slightly 
lower than 2,850 acres as estimated from the Landsat data. The difference seems to 
be negligible in practical application. Furthermore, Connor and Belanger (1981) 
reported that the Lake Washington water surface area was 2,851 acres at the stage 
14.6 ft. But, the surface area estimated from the Figure 5 is 2816 acres at 


the same stage of 14.6 ft. The difference Is about 1.2%. This difference Is 
hard to explain due to the lack of details associated with .the area computation 
method used by Connor and Belanger (1981). However, this 1.2% difference Is 
considered within limits of an acceptable deviation from the standpoint .of 
practical application. 

Lake Harris ; The dates were rearranged as shown In Table 2 In the In- 
creasing order of the lake stage. The water surface area at those four. dates 
was estimated from the Image 100 using the same procedure as used for Lake 
Washington. A typical color-coded classification result of water surface area 
for the Lake Harris on date 1/22/76 Is shown In Figure 6. The results of water 
surface area estimated by Landsat data on those four dates are shown In Table 2. 
The surface area varied from 17,430 acres at the. stage of 62.38 ft above the 
mean sea level to 18,657 acres at the stage 63.30 ft. The stages are plotted 
against the water surface area as shown in Figure 7. As Figure 7 shows, the 
range of stage fluctuation during the study period was only about 1 foot. This 
was mainly because the stages of Lake Harris are regulated well by the St. Johns 
River Water Management District. 

Storage Volume Determination 

Lake Washington ; The volume of lake was computed based on the method given 
in equation 1. The initial volume of lake at the stage 10.60 ft was not 
available. A constant volume designated as “C was used as a lake volume below 
the stage of 10.60 ft. The volume Increments at each stage of measurement were 
estimated. The results of the computed volumes are also shown in Table 1. 

The volume between the stages 10.6 and 15.85 ft was about 14,352 acre-ft. The 
lake volumes associated with each 0.01 ft lake stage Increment are also listed 
In Table 1. The lake volume changes ranged from 26 to 28 acre-ft for each 0.01 
ft increment In lake stage. The rate of change in Lake Washington volume in re- 
lation to lake stage Increments appears to be a stablized condition. This implies 


that the lake volume at other stages including stages below 10.60 ft, or above 
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15.85 ft could be roughly estimated based on the lake, volume stage Increment 
relation obitatned In this study. The computed lake volumes are also plotted 
against the lake stages as shown In Figure 8- This figure can be used to 
estimate the surface area associated with any specified lake stage between 

the measurements. 

Lake Harris ; The same technique used In Lake Washington was also used 
In Lake Harris, A constant volume designated as “C was used as a lake 
volume below the stage of 62.38 ft. The volume Increments at each stage of 
measurement were computed based on the method given In equation 1, The re- 
sults of the computed volume, and the lake volume associated with each 0.01 ft 
lake stage change are also listed In Table 2. The lake volume increment ranged 
from 176 to 183 acre-ft per 0.01 ft change In lake stage. This rate of change 
Is considered to be quite a stabllzed condition and can be used to estimate 
the lake volumes for other stages Including those below 62.38 ft, or above 63.30 
ft. The computed lake volumes were also plotted against the lake stages as shown 
In Figure 9 which can be used to roughly estimate the lake volume associated 
with any specified lake stage between the measurements. 

Further Applications 

Based on the conclusions made from previous sections, the technique 
developed in this study can be further applied In the following two cases. 

First, when the historical stage records are available, the historical 
Landsat data can be used to measure the water surface area associated with a 
certain stage recorded. Consequently, the relationship between lake surface 
area and lake stage can be established. The lake volume as related to the lake 
stage can also be estimated. Thus the district can directly correlate change 
in lake stage with the available water volume. Second, the technique can be 
applied to the case where historical stage records are not available, but the 
current stage record Is available or will be available In the future. The 


same technique as used In case- 1 can be used to estimate both lake water 
surface area and lake volume as related to the lake stage based on the 
current and future Landsat data and stage records. Thus, the historical 
lake stage and lake volume can be predicted based on the historical lake 
water surface area which was estimated from the historical Landsat data. 

Finally, an Important action needs to be taken by the district Is to 
establish a stage recorder for the lake that does not have any stage Informa- 
tion. Thus, the application In above case 2 can be accomplished after.. the 

stage data is recorded. 
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Table 1. Lake Washington volumo computation based on the water surface 
area measured from eight dates of Lanosat data. 


Date 

Water 
stage . 

Water 

surface area 

Lake Volume 
Different 5.01 ft 
stages Increment 


mm 

— Acre-*-** 


6/15/74 

10.60 

2.537 

C* 

26 

4/11/76 

11,81 

2.670 

C+3.150 

27 

2/9/76 

12.06 

2.683 

C+3,819 

27 

3/17/74 

12.45 

2.705 

C+4.870 

27 

2/27/74 

12.75 

2.729 

C+5.685 

28 

9/6/72 

14.80 

2.819 

C+U.372 

28 

11/29/73 

14.84 

2.827 

C+11.482 

28 

10/19/74 

15.85 

2.850 

C+14.352 



* C Is a constant volume below the stage of 10.60 ft 
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Table 2. Lake Harris volume computation based on the water surface • 
measured by four dates of Landsat data. 


Date 

Water 

stage 

Water 

surface area 

Lake Volume 

M'fferenr — OTTT 

stages Increment 


-Ft- 

—Acre— 


9/6/72 

62.38 

17.430 

c* 

176 

8/31/75 

62.50 

17.724 

C+2.109 

178 

2/14/75 

62.66 

17.963 

C+4.964 

183 

1/22/76 

63.30 

18.657 

C-t-16.682 



*C is a constant volume below -the stage of 62.38 ft. 



Figure 1. St. Johns River Water Management District 
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Figure 5. Relationship between water stage and water surface area for Lake Washington 
estimated from Landsat data. 
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Figure 6. Color-coded classification result of water surface area for Lake Harris, Florida 
on January 22, 1976. 
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I* Introduction to Landsat Remote Sensing System 

Landsats I, 2, and 3 wero launched by National Aeronautics and Space 
Administration (NASA) on July 23, 1972, January 22, 1975, and March 5, 

197S, respectively, to gather data about the earth's surface and telemeter 
those data to ground receiving stations. The operation of Landsat 1 was 
ended on January 6, 1978. The Landsat 3 follows- the Landsat 2 by 9 days. 
Landsat 4 will be launched sometime during the summer of 1982. 

Landsats 1, 2, and 3 were launched .Into a circular, near-polar orbit 
at an altitude of about 919 km. They circle the earth every 103.27 minutes, 
or roughly 14 times a day. Every 18 days the satellite passes over the 
same swath (185 km wide) of the ground. Landsat covers the earth's entire 
surface area except the northernmost 10% and southernmost 10% area. The . 
orbital Inclination angle of 99.11° was selected that provides a "sun- 
synchronous" orbit, so data collected are always at the same local time over 
the ssmie region. At higher northern latitudes the local time will be a few 
minutes late as compared to the equatorial crossing time. Meanwhile, local 
time at lower southern latitudes will be a few minutes ahead. There Is also 
a variation associated with the longitude, mainly because of the discrete 
time zones used. The average equatorial crossing time for Landsat 1 was 9:42 
a.m. till March 1977, then It was changed to 8:30 a.m. Landsat 2 was Initi- 
ated with an equatorial crossing time of about 9 a.m. and then was changed to 
9:30 a.m. In December 1977. 

Landsat 4 will be slightly different from the first three Landsats. For 
Instance, the orbital altitude Is 705 km and the earch coverage period Is 16 
days. 

Landsat program Is Intended to establish the usefulness of relatively 
reflective multi spectral Imagery for earth resources study. The application 
to the following areas has been Investigated: 




(1) agH cultural production; 

(2) water resource!? planning and management; 

(3) rangeland management; 

(4) forestry management; 

(5) land use and urban and regional planning and management; 

(6) environmental conservation and management; 

(7) geologic survey and mineral resource management; 

(8) marine resource, oceanography, and coastal engineering; 

(9) cartrographic and thematic mapping applications; and 
0.0) disaster warning and relief options. 
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I I . Sensor, Band Designation. Wavelength, and Resolution ; 

The typo of sensor, band designation, spectral sensitivity range, and 
spatial resolution used In Uandsat are shown In Table 1. Throe basic types 
of sensors are Involved: 

A. Return Beam Vidicon Camera : A multi spectral return-beam vtdtcon (RBV) 
television system, that operates In three spectral pass bands designated as 
band 1, 2, and 3 Is used on Landsats- 1 and 2. The wavelengths associated 
with the bands 1, 2, and 3 are 0.475-0.575 um (blue-green), 0.580-0.680 ym 
(yellow-red), and 0.690-0.830 ym (red-IR). The spatial resolution Is 76 m. 

Two panchromatic RBV cameras which replaced the three multispectral 
RBV cameras on Landsats 1 and 2, are being used on Landsat 3. The wavelength 
associated with this system ranges from 0.505 to 0.750 micrometer (yro), and 
the resolution is 40 m. This system has linear resolution which Is twice 
as better and still provides approximately the same ..areal, coverage as com- 
pared with the multispectral RBV cameras. 

B. Multispectral Scanner System : A Multispectral Scanner System (MSS), that 
operates in four spectral passbands designated as bands 4, 5, 6, and 7, is used 
in Landsats 1, 2, and 3. Furthermore, band 8 Is also used in Landsat 3 system. 

The wavelengths associated with the bands 4, 5, 6, 7 and 8 are 0. 5-0.6 ym 
(green), 0.6-0. 7 ym (red), 0.7-0. 8 ym (near IR), 0.8-1. 1 ym (near IR), and 10.4- 
12.6 ym (thermal IR), respectively. The resolutions for bands 4, 5, 6 and 7 are 
76 m and for band 8 is 234 m. At the time of this writing (March, 1982), the band 
8 Is not operating satisfactorily. . 

Typical representations of the four bands of the Landsat Multispectral 
Scanner (MSS) are used In this study. The characteristics of each band are 
listed as follows: 1 

Band 4: The band 4 emphasizes conditions of surface water bodies such as 

sediment loads, shallowness, directions and relative rates of flow. 




Band 5: The band 5 emphasizes cultural features- and exposed soil and rock 
surfaces and assists In the analysis of surface water conditions. 

Band 6: The band 6 helps to distinguish between live vegetation and land- 
water boundaries. 

Band 7: The band 7 penetrates best through haze and light clouds and also 
emphasizes live vegetation and land-water boundaries. 

C. Thematle Mapper ; A thematic mapper (TM), that operates In seven passbands 
designated as bands 1> 2 , 3, 4, 5, 6» and. 7 will be used on Landsat 4. The 
wavelength associated with, the bands 1, 2 , 3, 4» 5, 6, and 7 are 0.45-0.52 ym 
(blue), 0.52-0.60 ym (green), 0.63-0.69 ym (red), 0.76-0.90 ym (near Infrared), 
1.55-1.75 ym (Intennedlate Infrared), 10.4-12.6 ym (thermal Infrared), and 
2.08-2.35 ym (Intermediate Infrared), respectively. The resolution Is 30 m 
In all the bands except that the band 6 has a resolution of 120 m« Because of 
the resolution Improvement and wider range of spectrum, the Landsat 4 Is ex- 
pected to be a powerful tool for a Landsat remote sensing application. 
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III. Qround Receiving Stations 

The U. S. National Aeronautics and Space Administration (NASA) operates 
three Landsat tracking and ground receiving stations located at Greenbelt, 
Maryland; Goldstone, CallforJila; and Fairbands, Alaska. The ground. area re" 
ceptlon coverage for any receiving station Is about 28x10 km . In addition 
to the U. S. receiving stations, existing foreign ground stations Include two 
In Canada and one each In Brazil, Italy, and Iran. Besides, Landsat data 
outside the direct reception coverage areas of the existing ground stations 
can be acquired and recorded on the wide, band video tape recorders and then 
subsequently transmitted to a ground receiving station when the satellite 
passes by within the range. 

After the data have been received, either by real-time transmission or 
by tape recorder play-back, ground receiving stations record the data as wide 
band video recordings. These recordings are shipped to the NASA Image Processing 
Facility (IPF) at Goddard Space Flight Center (6SFC). Video tapes from the 
three U.S. stations are used to produce master 70 mm Images from all Landsat data 
collected by the multispectral scanner (MSS) for Landsatsl, 2, and 3. Second 
generation negatives of all Landsat scenes are furnished routinely to three U.S. 
Government distribution centers. Computer compatible tapes are produced by the 
IPF bnly when user requests are received at the USGS EROS Data Center. 


IV, General Electric's "Imago lOQ” System 

The General Electric “Imago 100“ Is a system which Incorporates both 
optical and electronic capabilities for enhancement purposes. This system 
Is one of the most advanced In photoanalyzers. The system Is designed to 
accept magnetic tapes with two tape drives and also contains a digitizer 
for use with films. A minicomputer associated with the system Is used to 
gather spectral data from known local areas and then compare data from other 
areas to decide whether the latter area meets the criteria as shown In the 
known, local area. Classifications appear graphically on the cathode ray 
tube display and outputs Includes the digital tape, colored CRT display, line 
printer, and color film recorder. 

Due to the very high price of this Instrument, only the limited facilities 
are available within the United States. These facilities Include such as USGS 
facility at Sixoux Falls, South Dakota; Jet Propulsion Laboratory In Pasadena, 
California; Johnson Manned Space Center, Huston, Texas; Goddard Space Flight 
Center facility In Beltsville, Maryland; and University of FlorWa-, Gainesville, 
Florida. 
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V. P. S, Distribution Centers 

Currently, Landsat Images are reproduced and distributed to users by 
three Federal Data Centers operated separately by the U. S'. Geological 
Survey (US6S), the U.S. Department of Agriculture (USDA), and the U.S. 
National Oceanic and Atmospheric Administration (NOAA). Tables 2 
and 3 list the U.S. Data Centers and the sources of available Image products. 
Image products provided by the data centers are similar In scale, format 
and type (Table 4). Interested users can request specific Information on 
characteristics of Image products available from the various data centers. 
USGS EROS Data Center does provide CCT's and various accession aids to users 
that are not available from the other. Federal data centers. 
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VI, The Northeast Regional Data Center. Florida 

The NEROC Is currently in the process of setting up procedures to handle 
digital imagery facilities for its user community. As one part of this, we 
are intersted in increasing the facilities and services available to users 
of Landsat data bases. 

LANDSAT DATA 

If you are a current or future user of Landsat machine-readable data from 
either the Multispectral Scanner (MSS) or Return Beam Vidicon (RBV) systems, 
NERDC would like you to share your data with us. 

Since the acquisition charge for the computer-compatible tapes is $200 
each, NERDC would like to make as many of these as possible available to all 
users. A copy of each tape that is made available will be placed in NERDC ?s 
tape library as a public reference volume available to all NERDC users. The ad' 
vantages of sharing this data are that: 

* It will make new data available to everyone; 

* It will provide you with a working copy of your data maintained by the 
NERDC at no cost to you. 
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V n . Uni versity of Florida, IFAS Remote Sensing Facilities 

The Uonornl Electric "Imaqe 100” located at NASA, Kennedy Snnce 
Center was relocated at the Institute of Food and Agricultural 
Sciences (IFAS) remote sensing system. University of Florida, In January, 
1982. About 1500 Landsat tapes came together with the Image 100 facility. 
This system will start operating sometime In the summer of 1982. The de- 
tailed Information on operation and personnel involved will be given to all 
potential users aruund- the state. 
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VIII. Landsat Newsletters 

"Landsat Data User NOTES*" a bi-monthly newsletter that presents 
Information about Landsat products, systems, and related remote sensing develop- 
ments, Is available to subscribers free of-charge.. Anyone, interested In 
receiving this publication should contact: 

User Services Section 

U.S. Geological Survey 

EROS Data Center 

Sioux Falls,. South Dakota 57198 

Another bi-monthly technical publication of Interest to Landsat users 
Is the "LANDSAT Newsletter," which provides both current statistics and historical 
data on satellite operations, as well as related Information on other earth survey 
programs. This newsletter Is available free of charge by contacting; 


LANDSAT Newsletter 
Missions Utilization Office 
NASA/GSFC 
Code 902.1. 

Greenbelt, Maryland 20771 
(301) 344-8826 

Users wishing to receive "Reflections," the quarterly newsletter of the 

Eastern Regional Remote Sensing Applications Center (ERRSAC) at NASA.'s 

Goddard Space Flight Center should write to: 

Reflections 
ERRSAC 
NASA/GSFC 
Code 902.1 

Greenbelt, Maryland 20771 
(800) 638-0748 

Users with Landsat data to share or users desiring Information about 
NERDC's Landsat computer support services should contact John Young at 
the NERDC or University of Florida, IFAS Remote Sensing Center. 
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IX. Landsat Update Information 

According to the recent (March. 1982) Information .update for us.ers of . 
National Oceanic and Atmospheric Administration's National Earth Satellite 
Service (NOAA/NESS) satellite data, several subjects as related to. Landsat 
are described as follows: 

Landsat 2 problem : The Landsat 2 yaw axis momentum wheel ha^ topped. 

This spacecraft has been the prime Landsat data collection platform. A pre- 
vious difficulty with this momentum wheel was overcome and It was restored 
to service in July, 1980. At the time of this writing (March, 1982.), the . 
satellite is not capable of collecting valid data. NASA is evaluating the 
problem and making efforts to restart the wheel, but without success so far. 
NASA regards the outlook for restoring the satellite to service as dim. 

Landsat 3 problem ; Variations in Landsat 3 multi spectral scanner mirror 
scan motions have been detected. They result in data loss on the left one- 
third of the picture. NASA has developed a ground processing correlation 
scheme and implemented it from the Goddard Space Flight Center facility to 
preserve the right two- thirds of the picture. Poleward of 35°N or 35°S, 
provides contiguous coverage. 

Landsat 4 system ; So far no problems or delays have come up that would 
interfere with the launch of Landsat 4 in the summer of 1982 or with the ini- 
tiation of NOAA's operational Landsat 4 program on January 31, 1983. 

The multi spectral scanner will be the operational sensor when NOAA first 
takes over the Landsat 4 system on January 31, 1983. It will take additional 
months to Implement routine Thematic Mapper operations. Presently, NOAA Is 
developing alternate Multispectral Scanner data acquisition schedules to find 
out which of them can be fitted to the engineering and cost constraints of the 
Landsat a system. 

System output will be capped at 136 processed Multispectral Scanner 


scenes per day for the first year or two of operation. Disaster and 
emergency events will have first call on this capacity and -soec1a1 acqui- 
sitions , for which requestors pay the full cost. Remaining capacity will 
be used to collect the Multlsoectral Scanner Basic Data Set ; scenes of 
general Interest acquired, to the extent possible, according to a published 
plan. The MSS Basic Data Set will deal with the routine scene collection 
objective of the operational system. 

Special Acquisition signifies Landsat 4 MSS scene data that are not 
scheduled for routine MSS Basic Data Set collection, but which are provided 
upon users' request according to the following five categories. 

1. Delivery of preprocessed digital data, to the requestor's site via 
communication satellite; provided MSS scene collected at a time and 
place specified by the requestor Is available. 

2. Delivery to the requestor of a frane of standard MSS imagery (not a 
color composite); If MSS scene collected at a time and place specified 
by the requestor Is available. 

3. Delivery to the requestor of a computer compatible or high density 
digital tape; If MSS scene collected at a time and place specified by 
the requestor Is available. 

4. Surcharge for delivery of a color composite to the user originally 
requesting the special acquisition of an MSS scene. 

5. Surcharge applied when the requestor establishes a maximum allowable 
cloud cover condition for the collection of an MSS scene. 


-13- 


Table 1. Campari sons of type of sensor, band designation , spectral 
sensitivity range, and spatial resolution used In different 
remote sensing systems 


Landsat 

Type of 
sensor 

Band (NASA 
deslanatlon) 

Spectral Spatial 

sensitivity range resolution 

1.2 

Return Beam 
Vidicon Camera 
(RBVl 

Band 1 
Band 2 
Band 3 

0.475-0.575 pm (blue-green 
0.580-0.680 wm (yellow-red) 
0.690-0.830 urn (red-IR) 

76 m . 
76 m 
76 m 


Multi spectral 
Scanner (MSS) 

Band 4 
Band 5 
Band 6 
Band 7 

0.5-0. 6 pm (green) 
0.6-0. 7 pm (red) 
0.7-0. 8 pm (near IR) 
0.8-1. 1 pm (near IR) 

76 m 
76 m 
76 m 
76 m 

3 

RBV 

two cameras 

0.505-0.750 urn (Panchromatic) 

40 m 


MSS 

Band 4 
Band 5 
Band 6 
Band 7 
Band 8 

0.5-0. 6 pm (green) 
0.6-0. 7 pm (red) 

0.7-0. 8 pm (near IR) 
0,8-1. 1 pm (near IR) 
10.4-12. 6um (thermal IR) 

76 m 
76 m 
76 ra 
76 m 
234 m 

4 

Thanatic 
Mapper (TM) 

Band 1 
Band 2 
Band 3 
Band 4 
Band 5 
Band 6 
Band 7 

0.45-0.52 pm (blue) 

0.52-0.60 pm (green) 

0.63-0.69 pm (red) 

0.76-0.90 pm (near IR) 
1.55-1.75 pm (intermediate IR) 
10.4-12.5 pm (thermal IR) 
2.08-2.35 urn (intermediate IRl 

30 m 
30 m 
30 m 
30 m 
30 m 
120 m 
30 m 


^Launch Dates: Landsat 1, 7/23/72 Operation ended 1/6/78 

Landsat 2, 1/22/75 
Landsat 3, 3/5/78 
Landsat 4, proposed Summer, 1982 


i 
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Table 2. Landsat data availability from U.S. Federal Data Centers. 



T 

U.S, Department of Agriculture 

Landsat Data 

USDA Aerial Photography 

Skylab Iraaoery 

\ 

■ 

National Oceanic and 

Landsat Data i 

; 

Atmospheric Administration 

NOAA Satellite Data 



Skylab Imaoery 


uses EROS Data Center Landsat Data 

NASA Research Photography 
USDA Aerial Mapping 


Photography 
Skylab, Apollo and 
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Table 3, EROS Data Center Imagery Holdings. 


Data Types 

Current Number 
of Frames 



Landsat I & II 
(MSS & RBV) 

1,067,000 



Sky lab, Apollo and 
Gemini Imagery 

50,800 



NASA Research 
Aircraft Imagery 

1,420,000 




USOI Aerial 
Mapping Photograph: 


3,123,000 
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Table 4. Sununafy of Landaat Image products available at EROS Data 
Center (ED). 


SIZE 

FILM PRODUCTS 

70 mm 

Black & white: positive and negative 
products (scale 1:3,369,000) 

9" X 9" 

Black & white: positive’ and negative 
products (scale 1:1,000,000) 

9" X 9" 

False color composite: positive 
product only (scale 1:1,000,000) 


PAPER PRODUCTS 

(Positive only) i 

’I 

9“ X 9" Black & white and false color composite . 

(1:1,000,000 scale) | 

20“ X 20" Black & white and false color composite ] 

(1:500,000 scale) 

i 

I 


40" X 40" 


Black & white and false color composite 
(1:250,000 -scale.) 
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ABSTRACT 

The methods of computing area from maps are classified under two . 
categories: (1) bounded area, and (2) weighted sub-areas. Pioneering 
methods of planiraeter* coordinate squares, triangular rule, trapezoidal 
rule, Simpson's rule, double-merldlan-dl stance (DMD), coordinates, and 
digitizer are compared with the relatively easy methods of weighing technique 
and finite segments, and also with a recently developed Monte Carlo 
process. 

The bounded area which Is used extensively In computing area of land 
use can be obtained by any one of three methods: graphical, arithmetical 
and computer. The graphical method Includes four techniques: weighing, 
planimeter, coordinates squares and digitizer. The accuracy of this graph- 
ical method Is highly dependent upon the skill of the analyst and the 
correctness of Instruments. The arithmetical method involves three rules: 
triangular, trapezoidal* and S1im)son. The accuracy of this arithmetical 
technique is dependent upon the number of offsets used to divide the entire 
region. If a region Is of very Irregular shape and of large size, then this 
arithmetical method Is very tedious and cumbersome to use. The techniques 
of double-merldlan-dl stance, coordinates, finite segments and Monte Carlo 
are Introduced by using a computerized procedure. The signs of latitudes, 
departures, and starting point are easily confused In the methods of OMD and 
coordinates. Finite segments can overcome some limitations such as a figure 
having some Inactive areas within the figure, but It Is also strongly dependent 
upon tl:e direction used to select the boundary nodes. Monte Carlo Is only 
a method which Is Independent of the direction and can be used In any compli- 


cated form of figure. 


The weighted sub-area which is used widely in computing the areas 
weighted to a measurement station can be obtained by any one of the four 
techniques: simple arithmetic methodi Thiessen method* stratified method 
and modified Monte Carlo method* however* the modified Manta Carlo method 
can be simulated directly by the computerized procedure and also has been 
applied successfully to compute the weighted area to each measured station. 

Finallyi from this study It can be concluded that the weighing technique 
Is relatively easy to use In the laboratory; the finite segments method Is 
a quick computerized procedure when the direction of selecting boundary- 
points are correctly performed; and the Monte Carlo process is more appllcabi 
and powerful than other computerized methods because this method Is not only 
Independent of the direction used to select the boundary nodes but also 
applicable to any complicated form of figures with arbitrary signs of lati- 
tudes* departures and starting points. 
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INTRODUCTION 

The computation of an area of ]and Is most frequently used 1n several 
fields of water resource systems. Area calculattons of Interest to water 
resource workers Include land use estimation, water quality, ecological 
systems, stratified sampling programs, etc. Pioneering methods used In 
computing area from maps, (hereafter tlie area computations are only re- 
ferred to the area computed from the maps), are those of planimeter, coordi- 
nate squares, triangular rule, trapezoidal rule, Simpson's rule, double- 
merldlan-dl stance (DMD), and coordinates (Drinker, 1969). Recently, a 
digitizer was also used (SAC, 1972). However, If an area Is very Irregular 
In shape and/or Is of large size or has Inactive areas within Its boundaries, 
then the previous methbds are very- difficult to adapt to a computerized 
procedure. Hence, there Is a need to Introduce some other practical methods 
from which the area can be calculated directly, either by the computerized 
procedure or by a weighing machine. 

OBJECTIVES ' ‘ “ 

The. purposes of this study are: 

(1) To discuss the existing methods used In computing areas; 

(2) To Introduce two relatively easy methods of weighing techniques 
and finite segments, and a recently developed Monte Carlo process, 
together with their ranges of application, and 

(3) To compare the advantages and disadvantages of each method. 

DESCRIPTION OF METHODS 

The methods of computing area from maps may be classified under two 
categories: (1) bounded area and (2) weighted to sub-areas. 
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BOUNDEO AREA 

The bounded, area has been widely used In determining the area of land 
Included within the boundaries. The area to be measured can be obtained by 
any one of the following methods: 

I. Graphical Methods: The area Is measured by using Instruments such as: 

planimeter, coordinate squares counter, digitizer and weighing machine. 

a. Weighing Method : A weighing machine used In a laboratory can be 
employed to obtain the area. Three procedures are Involved: First, 
a ratio between a known area and weight must be procured from a 
control unit. Second, the weight of an unknown area Is weighed. 
Third, multiplying the weight of an unknown area by the ratio which 
Is obtained In the first procedure gives the area of an unknown 
figure. The nonuniform thickness of the paper and changing humidity 
can affect significantly the accuracy of measuranent. 

b. Planimeter : The planimeter Is the commonest way of checking the area 
of a figure. It Is a small Instrument consisting of an arm, carrying 
a tracing point, which Is moved over the outline of the figure to be 
computed. Poor setting of the planimeter scale bar, and failure to 
check the scale constant by tracing a known area, can cause an error 
of measurement. 

c. Coordinate Souares : The figure Is marked off In squares of unit area 
The number of complete unit square Included Is counted, and the sum 
of the partial units are also estimated. A transparent paper marked 
in squares to some scale is placed over the figure and the number of 
souares and partial units counted. The number of squares can also be 
counted by a mechanical dot counter or a transparent paper dotter. 


Using coordinate squares which are too large makes It difficu 
estimte.the partial blocks and could cause an error in computation. 

d. molticer : A digitizer Is a device used to convert Information In 

graphic fom Into numerical Intelligence suitable tor processing, 
recording, or transmission on a digital data system. After the (x. y) 
coordinates of each point on the outline of the figure Is recorded, 
the size of the plane area also can be computed. This method Is highly 
dependent upon the direction used to digitize the boundary Coordinates. 

. Arithmetic Methods: A figure can be divided Into geometrical .shapes 

(triangles, trapezoids, and rectangles), and the follovring rules can be 

used to compute the area: 

a. Trlenoular Rule : A figure may be divided Into simple triangles.. The, . 

area can be computed by the f ormul a . 


Area » / S(S-a)(S-b) (S-cT 


where 

a, b, and c are the. sides of the triangle 

( 2 ) 

S = 0.5(a + b + c) 

Another formula Is used «hen an angle between two sides Is known. 

Area »0.5absinC 

where C Is the angle Included between two sides a and b. 
b. Tw.n.zo1dal Rule : If the figure 1s considered as made up of a series 

of trapezoids, all having the same base, the area can be determined 

based on the formula 


Area - | (h^ + 2 z h + hj,) 


where 


stance between offsets and 
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hg, h|j, and h ■ first, last and Intermediate offsets, 
c. Simpson *s Rule ; For generally parabolic areas, Simpson's one-third 
rale as follows Is applicable to obtain the size of the plane area. 

Ar«> • i (h, ♦ 2 £ ♦ 4 E ♦ h^) (6) 

where 

d • common distance between offsets, 

ha,hb,h^dd ^even " offsets. 

Furthermore, If the figure Is very Irregular In shape, then the 
arithmetical method Is very difficult to adapt by a computerized procedure. 
Therefore, the following computerized methods- are introduced; 

III. Computer Methods: The following four methods can be performed by the 
computerized procedure. 

a. Double-Merldlan-DI stance (DMD): The DMD of a traverse line Is twice 
the distance from a meridian through the most westerly station to the 
middle point of the line. The double areas of all of the trapezoids 
may now be found by multiplying each DMD by the adjusted latitude 

of that side. The area obtained by plus latitudes and minus latitudes* 
should be considered as a positive and negative area, respectively. 

The sum of these double areas will be double thd area of the figures to 
be measured. The disadvantages of this DMD method Is that the signs 
of OMD's, latitudes, departures, or areas are easily confused. 

b. Coordinates ; The area is equal to one-half the sum of the products 
obtained by each x-coordinate by the difference between the adjacent 
y-coordinate, taken In the same order around the figure. Similar to 
the DMD method', the signs of coordinates, latitudes, departures, or 
areas are easily confused. 
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If the signs and starting point are recorded carefully and 
without any Inactive areas within the figure, the above two methods 
can be used quite satisfactorily by using a computerized procedure; 

however. In practical application the above methods can accrue 

difficulty such as an Inactive area Included within the figure. 
Therefore, the following two methods which can overcome the 
limitations are introduced. 

c. . . Finite Segments ; The boundary of a figure Is defined by a series of 

linear segments between node points. The clockwise or counter- 
clockwise direction used to record (x, y) coordinates is dependent 
upon the figure either to be active or Inactive, respectively. The 
area is obtained by the following formula: 

M-1 y. + Y4 Yi + Ym 

Ar«i • (x,^, - X,) L+ (X, - x„) — ^ 

where X.j, Y.j * coordinates of boundary nodes; 

1 > boundary segment Index; and 
M = total number of segments. 

d. Monte Carlo Method: The Monte Carlo method is a procedure which takes 

advantage of the high speed of an electronic computer in solving complex 
problems In physical and mathanatical fields, htonte Carlo applications 
In the field of science and engineering are summarized In books by 
Hammersley and Handscombe (1964) and Shrleder (1967). This method 
Involves enclosing the figure to be measured within a rectangular 

area and by generating random numbers choosing points randomly distri- 
buted throughout this rectangle. The proportion of these points 
falling within the figured area Is, In the limit, the proportion of 
the rectangular area contained within the figure. But, in practical 
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application, the previous method has some limitations which were 
removed by Shih and Hamrick (1974). 


^ • To Determine Whether a Random Point Is Within or Without the 
watershed! 


Shih and Hamrick (1974) developed an alternative test 
based on the principle that for any completely bounded region, 
a radial line constructed in any direction from a given point 
must cross the boundary an odd number of times if the point Is 
located within the region or an even number of times 


If without (assuming zero to be an even number). This test Is 
ambiguous only In the case of node points (the intersections 
of straight line segments representing the figure boundary). 

They also developed a second rule with a computerized technique 
for solving the ambiguity that exists when the radial line 
penetrates the boundary at the node point. 

11 . Procedure of Computation ; 

Based on the technique described in the previous section, 
the following procedures are used to obtain an analog of any 
shape of figure. 

(1) Enclose the irregularly shaped boundary with a rectangle 
whose coordinates are also recorded. 

(2) Read the X and Y coordinates of the boundary segments. 

(3) Compute the weighting factor of each boundary node according 
to Rule 2 (Shih and Hamrick, 1974). 

(4) Generate random points with uniform probability over the 
enclosing rectanole. 


-9- 


(5) Draw an Imaginary Une from the random point and parallel 
to the x-ax1s. 

(6) Count the number of intersections of this line with the 
boundary 

(7) Test whether the random point falls within the boundary 
according to Rule 1 (an odd or even number of inter- 
sections). 

(8) If the above test succeeds, increase the counter of accepted 
points by one. 

(9) Repeat the processes of 5, 6, 7 and 8 until the number of points 
assigned is reached. 

(10) Compute the area by dividing the accepted points by the total 
number of random points, and multiplying this ratio by 
the enclosing rectangular area. 

The above procedure is simulated as a flow chart shown in Figure 1. 
iii . Selecting Boundary Segments 

The boundaries of a watershed can be defined by the given 
coordinates of successive points along the boundary (in a clockwise 
direction) and considering the boundary between each pair of successive 
points to be a straight-linear segment. The actual boundaries could 
bo approximated as closely as desired by increasing the number of 
such- segments ; but, the user should note that the more segments chosen 
the more computing time and user's time are required. A later 
example will show the effects of many irregularities of natural 
boundaries such as lake shorelines and the watershed may be averaged 
by reducing the boundary planform to a simplified polygon. The 
general principle of this procedure is to represent the planform 
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Figure 1; Flowchart for Solving Integrated Area 




I 

.} 

with as few sides as possible without changing the basic shape of 
the boundary. 

1v.. Convergence of Weights 

Because the Monte Carlo method relies on the laws of probability, 
a large number of random trials should be taken In the Interest 
of precision. The method normally used to estimate the relationship 
of sample size and accuracy Is the large-sample normal approxi- 
mation. Using the Central -Limit Theorem, the binomial distribution 
can be approximated by a normal distribution for a large N, where 
N Is the total number of trials. The sampling error of any statistic 
Is proportional to l/(M)^/2^ Tl,e convergence of the Thiessen 
weights Is a statistical convergence I.e., the probable error of 
estimation Is proportional to 1/(N)^^^. A detailed discussion was- 
given by Shih and Hamrick (1974). 

WEIGHTED SUB-AREAS 

The average amounts of environn«ntal elements such as water quality, 
ecology system, and land use, etc. over a specific area is required In many 
water resource problems. Thus, mean value problems can be solved by sampling 
techniques. For convenience let variables A-j , . . . , Aj^, and , . . . , X|^ be the 
subareas and measured values of stations 1, ...k, respectively. Then the 
estimated weighted-average amounts for the region are 

X'WlXi + ••• (7) 

where « A^/A, weighted area; and 
A ■ A-j + ... ♦ A|^, total area. 

In most cases the values of X^ are obtained first by laboratory experiments 
or field measurements, and then the values of are estimated based on the 
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following four methods: 

I. Simple Arithmetical Method: The values of are assumed equal to one. 
Equation 6 can be rewritten as 

X • z X./k (8) 

1»1 ’ 

Equation 8 Is the simplest method which can give a good estimation of 
average value In a flat area under the condition of measurements that 
are uniformly distributed and the Individual measurement does not vary 
widely from the true average. However, this method does not take into 
account the measurements outside, but near the boundaries of the area. 

II. Thiessen Polygon Method: Thiessen (1911) developed a method which attempts 
to allow for nonuniform distribution of measurements by providing a 
weighting factor for each measuranent. The measured points are plotted 

on a map, and connecting lines are drawn. Perpendicular bisectors of 
these connecting lines form polygons around each measured point. *nie 
sides of each polygon are the boundaries of the effective area assumed 
for the measured point. The values of Aj are determined by the methods 
indicated In the section of bounded area computation. However, the 
limitation of this method Is Its inflexibility; for Instance, a new 
polygon, being reoulred every time there Is a change In sampling location. 
Also, the method makes no attempt of overcome the orographic Influences. 

III. Stratified Method: This method Is a plan by which the region Is divided 
Into homogeneous subregions or strata. In computing a strata the analyst 
can make full use of his knowledge of orographic effects or other 
Influential factors. After constructing the strata, the following equation 
Is used to calculate the mean values, X, 
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where 


M 

2 

1»1 
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( 9 ) 


M ■ total number of strata 

■ total number of observed stations In the 1th stratum 
1,j » Index of stratum and observed station, respectively; and 
» observed value of the jth station In the 1th stratum 
The values of can be determined by either simple arithmetical method 
or the Thiessan method. If a simple arithmetical method Is used, then 
equation 9 Is reduced as 

X- I VI.X^ tio) 

i“l 


where 

» A^/A, weighted area of the 1th stratum; and 


X, » 2^ Xii/kj, average observed area of 1th stratum. 

1 j»l ^ 

If a Thiessen method is used to perform the subpolygon for each station 
then equation 9 is used to calculate the weighted area. The greatest 
limitation of the stratified method is also Its inflexibility. 

IV. Modified Monte Carlo Method: As Indicated in previous sections, the 
Thiessen polygon and stratified methods suffer from their Inflexibility 
in that a new Thiessen diagram is required every time there is a change 
in the sampling location or a recorded station with missing data. This 
limitation can be overcome by the modified Monte Carlo method (Shih and 
Hamrick, 1975). 
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a . Proceduros of Computation ■ ' 

The following procedures are used to compute an area weighted to 

each measured station: 

1. Compute the weighting factor of each boundary node which relies 
on the procedures 1, 2 and 3 as indicated In the jrrevlous section 
of. the Monte Carlo method. 

2. Determine whether a. random point falls within any shape boundary 

which Is based on the procedures 4» 5» 6, 7 and 8 as Indicated 

In the previous section of Monte Carlo method.. | 

3. Assign the undora point which Is falling within the boundary to 

the nearest measuring point. 

4. Repeat processes 2 and 3 until a predetermined large number of 
points are reached. 

5. Compute the relative area ratio of the bounded region to the 
enclosing rectangle by dividing the number of accepted points by 
the total number of random points. 

6. Calculate the computed weights of each measuring point by dividing 
the number of points assigned to each measuring point by the 
total number of accepted points. 

7. Check whether the sampling location Is changing^ 

8. If the response to 7 Is yes. the processes from 2 to 7 are repeated. 

The above procedures are simulated In. a flow . chart In Figure 2. . 

b. An Irregular or L-Shaoed Watershed 

An expected accuracy of computation by this Monte Carlo method 
depends upon not only the number J random points, but also the shape ( 

of the watershed. The number of random points that affects the sample 
error has been discussed In a previous section. However, the 
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efficiency of this Monte Carlo method affected by the shape formed 
by the watershed. should also bo discussed because, a large number of 
random points fall off the outside boundaries of the watershed 
which differs greatly from that of the enclosing rectangle. This 
difficulty can.be avercome..by following three techniques; 

i, Eoual Rectangles ; A watershed Is enclosed by a number .of 
smaller rectangles of equal area that have, a common edge 
which cuts the watershed. In order for this method to be 
used more widely* the following relationship should be 
Introduced. Let A.|, ..., Aj, be the relative area ratio 
falling within the boundaries of sub-watershed 1, ...» n. 
Then the new relative area ratio, R] , .... R^» of sub- 
watershed 1 , ...» n are equal to 


An/ i A. 
‘ 1=1 ^ 


^ • V L/i 


1=1 


(11) 


The final computed weights, Wi , ..., W^, of rainfall stations 
1 , ... , m are equal to 

. . 

: 02) 


n 


Wm • E R.E. 

"• 1 1m 


Where E.jj Is the computed weight of rainfall station j In 
subrectangle 1; j Includes the rainfall station from 1 to 
m and 1 includes the subrectangle from 1 to n. For example 
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11 . 


111 . 


m Is the total ..Ruinber of rainfall stations and n Is the total 
number of subrectangles. 

Unequal Rectangles ; The technique of unequal- rectangles Is 
similar to the equal rectangles method except that the water- 
shed Is. enclosed by a number of. smaller unequal rectangles. 

Let S.|, •••♦ represent the area of enclosing rectangles 
1, ..., n; Ai, ...A^ and R-|, ...,R„ are defined In the case 
of equal rectangles. The value of R.j Is 
n 


Rl « A^S^/ E A^S^ 

« I ♦ 

• • 

• n • 


(13) 


1® I 

The final computed weights, W<| , ..., of rainfall stations, 

1, ..., m, are similar to equation 12, except that the R.j 
values are replaced by equation 13. 

Single Rectangle ; The more random points chosen, the greater 
the accuracy of the estimates obtained. Therefore, in a 
watershed which has a lower relative area ratio, the single- 
rectangle technique Is still applicable by Increasing the 
random trials. A detailed technique of application was given 


by Shih and Hamrick (1975). 


Hew Thiessen Coefficients for Missing Data 

As Linsley et al. (1958) Indicated the greatest limitation of the 
Thiessen method is Its Inflexibility, because a new Thiessen polygon 
is required every time there is a change in the gage network. This 
modified Monte Carlo method can be used to overcome this limitation. 


In general, there are two cases of missing data. Case 1: The missing 
data of each rainfall station are priorly known, and any missing periods 
of record are assigned as a new station set. The distance of a random 
point from all rain measuring stations Is calculated simultaneously In 
each station set, and the random point Is assigned to the nearest rain 
measuring station In each set. Case 2: The missing data of each rain- 
fall station Is posteriorly known, I.e., how many stations with missing 
data are unknown. In this case. If a station with a missing record Is 
found, then that station Is omitted and a new gage network Is considered. 
Based on this new gage network, a computed weight Is performed by a re- 
peating procedure. A detailed description of these procedures will be 
discussed In the section of computer program. 

COMPUTER PROGRAM 

Based on Equation 1 of finite segments, procedures for computation In 
Monte Carlo and modified Monte Carlo methods, a systematic flow chart for 
the computer program development Is shown In Annex 1. Nomenclatures for 
the computer program are listed In Annex 2. The systematic flow chart 
was converted to a computer program for the CDC 3100 computer with Fortran 
IV language. The users manual for the COC 3100 program Is also presented 
In Annex 3. The computer programs are listed In Annex 4. If the other 
computer systems are used, some -mod4.fi cations of the computer program 
might be needed. 


EXAMPLE OF APPLICATION 

A example Is used to Illustrate the application of computing the 
bounded and weighted sub-areas. 
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In order to test the applicability of these newly developod techniques, 
an irregular area as shown in. .Figure 3 not only has inactive areas within the 
figure, but also can be calculated easily by the triangular rule. As mentioned 
in the previous sections, the direction used in setting boundary nodes is a 
“very important feature in some methods, because an improper procedure can cause 
a serious error. In order to investigate this nature, as Figure 3 illustrates., 
four different sequences used to select the coordinates are demonstrated: 

Case lA: 1-2-3-4-5-6-7-8-9-10-6-11-12-13-14-11-5-1. 

Case 2A: 1-2-3-4-5-6-10-9-8-7-6-11-12-13-14-11-5-1, 

Case IB: 1-2-3-4-5-6-7-8-9-10-11-12-13-10-9-14-6-5-1, 

Case 2B: 1-2-3-4-5-6-14-9-10-13-12-11-10-9-8-7-6-5-1. 

Based on these four cases, several methods are used to obtain the areas. The 

results are shown in Table 1. As can be seen from Table 1, the cases of lA 
and IB have a serious error by using the methods of finite segments and 
digitizer. The results of areas calculated by the Monte Carlo method with 
2000 random walks are all the same values in four cases. This implies that 
the technique used to select the consecutive points along the boundary in either 
a clockwise or counterclockwise direction, or a combination of both directions, 
gives no difference in thr results obtained by using the Monte Carlo method. 

This is a very useful tool when the directions used to select the boundary 
nodes are mixed up in both directions. The accuracy of Simpson's method is 
directly proportional to the number of divisions. For example, the number of 
eight and ten offsets as shown in Figure 3 are called Cases 1 and 2, respectively. 
The calculating results are also shown in Table 1. As Table 1 shows, Case 2 
has a better solution than Case 1. The results obtained by the coordinates 
and weighed method give a good agreement with the results obtained by the 
triangular rule. 
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F1g. 3. Area Under an Irregularly Shaped Form. 





Diff. 

Cases 

Graphical Method 


Arithmetical 

Method 

Computerized Method 

Coordinate Digitizer 
Squares Method 

Wei ghed 
flethod 

Triangular 

Rule 

Simp. 

Rule 

Monte 

Carlo 

Finite 

Segments 

lA 

16.97 34.03 

17.18 

17 

18.3 

17.06 

35 


2A 16.97 17.03 17.18 17 17 17.06 17 


IB 20.84 29.82 19.98 21 20 20.81 29 
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RESULTS AND DISCUSSION 

After applying the techniques to several practical problems, the following 
con clxis Ions were drawn: 

1. The accuracy of the graphical method is highly dependent upon the 

skill of the analyst and the correctness of the instruments. Although 
the results obtained by the coordinate squares, pi ani meter and 
weighing technique are In good agreement with the results obtained 
by the triangular rule, some observations must be mentioned. For 
example, the nonuniform thickness of the paper and varying humidity 
can significantly affect the accuracy of weighing measurement. Poor 
setting of the planimeter scale bar, and failure to check for the 
scale constant by tracing a known area. c?.n cause an error of 
planimeter measurement. Using coordinates which are too large make 
it difficult to estimate the partial blocks from which an error of 
computation can occur. The cases of lA and IB will have a serious 
error if the digitizer is used. These results show that the digitizer 
Is highly dependent upon the direction used to digitize the boundary 
coordi nates . 

2. The accuracy of the arithmetical technique is dependent upon the number 
of offsets used to divide the entire region. Case 2 with ten offsets 
has a better solution than Case 1 with only eight offsets. If a 
region has a very irregular shape and large size, then this arith- 
metical method is very tedious and cumbersome to use. 

3. The signs of latitudes, departures, and starting points are easily 
confused in the methods of DMD and coordi natef. Although finite 
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segments can overcome some limitations such as a figure having some 
Inactive areas within the figure, as Case lA and IB shows; It Is strongly 
dependent upon the direction used to select the boundary nodes but 
also applicable to any complicated form of figures with arbitrary sign 
of latitudes, departures and starting points. 

4. The modified Monte Carlo method has been successfully applied to compute 
the weighted sub-area. The result also indicated that only the modified 
Monte Carlo method can be simulated directly by a computerized procedure 
to compute the weighted sub-area to each measured station. 

5. The modified Monte Carlo method can be applied to compute the weighted 
sub-area not only for rainfall stations but also for any type of 
estimating mean value over a specific area. For example, an average 
of nutrient concentration within the Take can also be computed based 
on the modified Monte Carlo method. 

6. The-‘modif<ied Monte Carlo method has been extended to compute the new 
weighting factors when data are missing. 
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CALL SUBROUTINE —121 
FINITE SEGMENTS 
METHOD USED TO 
COMPUTE AREA 





Annex 2 . Nomenclature For Computer Programs. 


Variables of Input: 

IN - Methodology Indicator, the finite segments method used 

1 ■ Finite Segments Method Used 

2 > Monte Carlo Method Used 

3 ■ Modified Monte Carlo Method Used 
0 ■ End of Job 

TITLE - Input data Identification 
N - The number of boundary points chosen 
M - The number of measuring stations. 

NSET - The number of Monte Carlo Points assigned 

XMIN - Mini ml m range In X axis 

XI4AX - Maximum range In X axis 

YMIN - Minimum range In Y axis 

YMAX - Maximum range In Y axis 

FAC - The scale of feet used In per unit of length, zero means 
the dimensionless of area 
X - X coordinate of the boundary point 

Y - Y coordinate of the boundary point 

AX - X coordinate of the measuring station 

AY - Y coordinate of the measuring station 

Variables of Output: 

AF - Relative area ratio 

AREA - Total area In acres 

I - Measuring station Identification 

WF(I) ♦ Thiessen coefficient of station I 

FACT - Area converting factor in unit of acres 

SUM - Total area in acres computed by finite segments method. 


Annex..3* — Users Manual 
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For Program to Compute the Areas and Thiessen 
Coefficients Based on Finite Segments, Monte 
Carlo and Modified Monte Carlo Methods 
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Program Limitations 

1. Limit of 400 boundary points por area 

2. Limit of 500 measuring station per area 

Reouisitlon for .Computer WOrk 

Estimated Time - 5 seconds Is needed In Finite Segment per area, 

15 seconds ^s needed In Monte Carlo method per 1000 random points. 

20 seconds Is needed In ^tod1f1ed Monte Carlo Method per 
1000 random points 

Category - Production run 


Disk - 6000 

FORMAT INFORMATION 

Symbols used to Indicate the proper method for numbers or letters entered 
In card columns shown are : 

RJ - Indicates that a whole Integel number must be right justified 
In- card columns shown 

DP - Indicates that- the number must have a decimal point Indicated in 
one of the card columns. 

A .-.-any alpha-numeric character. 


CARD FORMAT INFORMATION 


First Card: Control and Parameters Card 


C.C. Symbol 

1-5 RJ 


A 


Description 

An Integral value Is required. 1, 2, and 
3 Indicate that the finite segments, Monte 
Carlo, and Modified Monte Carlo are used, 
respectively. 

Title Input data Identification. 
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Second Cav'd : Parameters Cav’d 
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c.c. 

Symbol 

Description 

1-10 

lype I: 
RJ 

when the integer 1 is shown in first card on 
column 5 

Total number of -boundary points 

11-20 

DP 

Scale of feet used in per unit of length 

1-10 

Type II: 
. M 

When the integer 2 is shown in first card on 

column 5 , ^ ^ a 

Total number of random points assigned 

11-20 

RJ 

Total number of random points assigned 

21-30 

DP 

fiinimum value of enclosing, rectangle 
coordinates along X-axis 

31-40 

DP 

Maximum value of enclosing rectangle 
coordinates along X-axis 

41-50 . 

.DP. 

Minimum value of enclosing recwngular 
coordinate along Y-axis 

51-60 

DP 

Maximum value of enclosing rectangular 
coordinate along Y-axis 

61-70 

DP 

Scale 4 feet used in per unit of length 

1-10 

Type III: When the integer 3 is shown in first card on 
column 5 . ^ 

RJ Total number of boundary points 

11-20 

RJ 

Total number of rain measuring stations. 
If a relative area ratio of study area 
to enclosing rectangle is expected, only 
the value of 1 should be used. 

21-30 

RJ 

Total number of random points assigned. 

31-40 

DP 

Minimum value of enclosing rectangular 
coordinate along X-axis 

41-50 

DP 

Maximum value of enclosing rectangular 
coordinate along x-axis 

51-60 

DP 

Minimum value of enclosing rectangular 
coordinate along y-axis 

61-70 

DP 

Maximum value of enclosing rectangular 
coordinate along Y-axis 

71 -80 

DP 

Scale of feet used in per unit of length 




Third Card; Boundary Points Coordinates Card(s) 


c.c. 

Symbol 

Description 

1-5 

DP 


11-15 

DP 


31-35 

DP 

X coordinate value of boundary segments 

41-45 

DP 

choosing in clockwise direction 

51-55 

DP 


61-65 

DP 


71-75 

DP 

• 

6-10 

DP 


16-20 

DP 


26-30 

DP 


36-40 

DP 


46-50 

DP 

y coordinate value of boundary segments 
choosing in clockwise direction. 

56-60 

DP 

.... 

66-70 

DP 


76-80 

DP 



Note: The maximum boundary points included per card is onl^y f 

the card can be used as much as required in number of stations, 
example, 7 cards are needed in 50 boundary points. 


Therefore, 

For 
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MeasuHna Stations Coordinate Card(a): This card Is required only when the 

^ Integer 3 Is shown In the first card 

on column .5 


c.c. 

Symbol 

Descriotlon 

1-5 

OP 


11-15 

DP 


21-25 

DP 


31-35 

DP 


41-45 

DP 

X coordinate of rain measuring station 

51-55 

DP 


61-65 

DP 


71-75 

DP 


6-10 

DP 


16-20 

DP 


26-30 

DP 


36-40 

DP 

^ y coordinate of rainfall measuring station 

46-50 

DP 

56-60 

DP 


66-70 

DP 


76-80 

DP 


Note 1: 

The maximum measuring stations Involved per card Is only 8 stations. 


Therefore, 

use as many new cards as necessary. 

Note 2: 

Use as many 

new control and title cards with succeeding cards as 

necessary. 

The last card must be present as a blank form. 


fW 




ANNEX 4. Computer Program with Fortran IV Language 
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ORIGINAL PAGE 18 
OF POOR QUALITY 


X IS TmE X COO><DINATE OF Tm£ BOUNDARY 
Y I« THE Y COO»OInaT£ OF IMF BOUNDARY 
NN IS A tdEIOHTiNG FACTOR Fo« NOOt POINT OF MOUNOAPY 

lA and IX APE Initial ood integer for using in ranou suRPouriNt 

M IS THE NUMBER OF BOUND ftPY POINTS CHOSEN 
M IS THE NUMBER OF RAIN mcasURING STATION 
NSET IS THF number OF RANDOM POINTS FXPErTEO 
XMIN and XMAX are THE m|nImUM AND HAkImIJm RANGF IN X AxIS 
YMIN AND YMAX AR£ THF MINIMUM ANO MAXIMUM RANGE IN Y AXIS 
AX IS the X COORDINATE OF Th£ RAIN mfASIJPING STATION 
AY IS THE Y COORDINATE OF THE RAIN MFAStIPiNO STATION 
FAC IS THF FEET USFO IN oi>’R UNIT 2£«0 mfaNS OIMENSIONLESS OF AP£a 
BFAL L(500) 

dimension X(AU0 ♦Y(hO*') •AXIS.'.'P) ♦AY(5C0-) ♦WF(5t)w.) ♦NSCSOC) •YY(400) ♦ 

- INN(AOO) •TITLF(O) 

NRs*S.) 

N-i»a6l 

READ INPUT DATA 
33 REAO(NR*3) IW^TITLE 
3 FOPMArdS.3x.MAo) 

IA*S 

ixs7 

CHECK vKHFThcR the FND OF «;TaTION SET 
IFdW.EO.c) GO TO rtS 
WRlTctNW.^) TITLE 
5 FORMAT (lril*//5X*MAB//) 

IFdW.EO.l > GO TO 1 IM 
iFdW.rO,?) GO TO d '' 

READ (NR to) NtMtNSET.XMlWfXMAXtYMlN.YMAX tFAC 
A FOPMAT(3n0.5FU.3) 

I F AC*F AC» 0 .01 
FACTsFAC*FAC/A35ft ).'* 

IF(IFAC.EO.n) FACTal.:- 

WRITE (NW. I an) Nt.MtNSET tXMlNt XMAX* YMlNt YMAX tFACT 
10« FORMAT (//) 3X.2MHNUMBER OF BOUNDARY SEGMENTS *. I4/1 0X.2'^MNUMB£R OF 
IWAINFALL STATIONS stI4/i Xt25MNUMB£R OF RANOOM POINTS stI6/lOXtl6H 
aX-AXIS minimum s.Frt.4/1 XtlbHX-AXlS MAXIMUM a t Ffi ,4/ 1 ») X ♦ IBHY-AX IS M 
3INIMUM 8.rB,4/lOXtloMY-AXlS MAXIMUM a , FB .4/ 1 C X .24HARF A . CONVERT ING 
4FACT.)P a. FIS. 5//) 

REAO(NRtM) (Xd ) tY d) t laitU) 

M format (INFS. 1) 

R£AD(NRt4) (AXd ) t AY d ) . laj tM) 

*» format (16FS. I ) 

AOn A SMALL value to EACH NODE FOR ORTA I MG THE NN FACTOR 
00 l?S 1*1 tN 
1D6 YY(I)aYd)tO.(;JI)''l 

COMPUTE THE WEIGHTING FACTOR NN FQR FACH BOUNDARY NODE 
DO r?l laj.N 
NN ( I ) *0 

IF(Yd).FO.Yd.l) » GO TD \ i,? 
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C 


103 

lv'2 

ni 

in? 


IV tI ^ { .j ‘ ^ .y ( n ) mnm ) snn (I ) m 
GO ^To •ANo.yy u ) . gt.y < n > nn( nswN( n ♦! 

NN < I ) « 1 
CONTINUE 
no n? 1*1 *M 

NS (I) = 0 
CONTINUE 

NAs(> 


NW*Vs 3 

XMNsXMAX-XMIN 

YMNsYMAX-YMIN 

no 500 iksunset 

IKSU 

IU*0 

GENEPATE THE RANDOM NUMrt£R 
CALL RAN00(IX»IY»R0M) 
iXsiY 


C 


XTaXM-lN ♦ «UM*XMM 
CALL RANOD ( IA« IH«»DN) 

IA*IB 

YTsYMIN ♦ R0N<»YMN 
X(N^l)sX(l) 

'■'■te^section ALOMG the X AXIS IN either side 
IF(YT.EO.Y(K).AN n,XT.EO.X(K)) 60 TO 31“. 

GO TO 10 


IC 


11 

U 

2( 

‘♦I 

<♦2 

30U 


IF(YT*EQ*Y (K) J GO TO ?-> 

IF(Y(K) .GT.YT.ANO.YT.GT.yiK^lJ ) 
IF(YT#GT#Y(K) ♦AND* Y (K^l ) «0T*Y**) 
GO To 300 

IF(X(K),LT»XT •AN0*X (K^ 1) ,GT«XT) 
IF(X(K) •GT.Xr.ANO*X(K^l ) ,LT,XT) 
IFUT-X(K)) IU.31C.,12 
IR*IR^NN(k) 

60 TO 3ft 0 
IL*IL«>NN(k) 

GO TO 30.“) 

IF (X (K)-XT) 


GO TO 40 
GO TO 40 


Go 

GO 


To no 
TO 310 


12v310«ll 

IR*I»*1 
60 TO 3(vft 

IL*IL»1 

continue 


31- nasna.i ^ ‘ 30E 
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SAVeaL(l) 

ISUH *1 

nn 91 I s?*M 

L(I)*(XT-AX(n )♦•<? ♦(YT-AY(I))**2 
IP‘(L(I)-'5AVE) 13. 91*91 
13 SAVFsL(I) 

ISUHsI 
91 CONTINUE 

NS(ISUR) *NS(ISUB)^l 
NO TO 

3:)2 KjppsNPW^l 
S'l CONTINUE 

C COMPUTE RELATIVE A»EA RATIO AWO <iEIOHTUK' FACTOR OF EACH STATION 
AF* IMA/ (Float (NA*NPR)) 

AREAaAF*(XMAX-XMIN)*<YMAX-YMIN) *F aCT 
C PRINT THE results 

wRITE(Nw.7> AP.APFA 

7 FORMAT (IOX.^IHRElATIVF AREA RATIO a*F9.6.2AH TOTAL AREA IN ACRES 
I a.FlS.6//) 

DO Z\ Isl*M 

<<F(I) sNS(I)/FlOAT(NA) 

^\ w»ITE(N*<.A) I.WF(I) 

A format ( nx.SSHCO'^PUTEO >*tIOHT OF RAInFALL STATION. I9.3H s.F9,6/) 

GO TO ft 3 

ID9 CALL ARAA(X.Y.NR.NW) .... 

GO TO A3 

Uv* call MONTF(X.Y.YY.NN.NR.MW) 

GO TO 83 
AS CALL EXIT 
FUO 


SUHROUTINE ARAA(X. Y.fJft.MW) 

DIMENSION Xin.YU) 

C RF.AO INPUT DATA 

RFAD(NR.A) N«FAC 

6 FORMAT(IIO.F13.3) 

REAO(NR.i»I (X(I) .Y(I> .Ul.N) 

4 FORMAT (IbFS. I ) 

IFACsFAC^C.Cl 
FACT=FAC*FAC/435A»'.r 
IFdFAC.CO.O FACT«1,{ 

^RITEINW.IOA) N.FACT 

1?A FORMAT (//I 3X.?9HNUMHER OF BOUNDARY SEGMENTS s, Ii./IOX.EAMAREA CuNVE 
IRTING FACTOR a.FlS.S/ZI 
SUMsii.O 
X (N*l > sX U ) 

Y (N*Hmy n ) 

no 7 Ksl.N 

SUMla(X(K.l)-A(K) ) ♦ ( Y ( K ♦ I ) *Y ( K )) /? . 0 

7 SUMsSUM*S,|MI 
SUMsSUM«FACT 

woiteinr.a) sum 
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“ FOf<MAT(/?X»aiHT0rAL APEA in acres a*FlH.6//) 

PFTUPW 

FNO 



SUHROUTINir monte (XiY*YY*NN»f4P»Nnn 

dimension X(l)*Y(ntYYn)«NN(l) 

TA*5 

IX*7 

C READ INPUT DATA 

PFAO(NW*H) N*NS£T«XMIN»xMAX*YMIN*YMAa*FAC 
d FORMAT (SIlOtSFlO.l) 

IFACaFAC^O.Ol 
FACTaFAC*FAC/435^0.C 
IF(IFAC*EO.(n FACTsl,<' 

WRITECNWflOH) N» NSET»XMIN»XMAX«YMIN,YMAX*FACT 
10ft FORMAT (//10X*?NHNUMHEP OF BOUNDARY SEGMENTS st I4/10X.25HNUHHER OF 
IRANOOM POINTS «♦ 16/lOX, 16HX-AAIS MlNlMl>M * ♦F8.4/1 OX* 16HX-AXI S rtAXi 
2MUM a.F8./*/lOX*l6HY-AXIS MINIMUM a ,F8. A/I OX ♦ 16HY-AX IS MAXIMUM a^pB 
3.4/nx*24HAREA converting FACTOR a,FlS.5//> 

REA0(NR*9) {X{r>«Y(n*Ial*N) 

W FOPMAT(lftFS.l) 

C ADD A small value TO EACH NODE FOR OMTAImG THE NN FACTOR 
DO 136 1»1*N 
IC.6 YYm*Y(I)^0.900»'l 

C COMPUTE THE WEIGHTING FACTOR NN FOB EACH BOUNDARY NODE 
DO lOI Ial*N 
NN<I)a9 

IF(Y(I) .EO*Y(I^I) ) GO TO 102 

IF(YYU) .LT.YU»l).AND.YY(n.GT.Y(I) ) NN(I)aNN(I)^l 
IF(I.EO.I) GO TO 103 

IF(YY(I) .LT.Yd-n .ANO.YYU) .GT.Yttn NNUlaNNd) *»I 
GO TO lOl 

i03 IF(YY(I).LT.Y(N).ANlJ.YYd).GT.Y<In NN ( I ) aNN (I » ♦!. 

GO TO 131 
H2 NNd)sl 
1 01 CONTINUE 

C GENERATE THE RANDOM POINT 
NAaO 
NWR»0 

XMNaXMAX-XMlN 

YMNaYMAX-YMIN 

DU SOP lK8l*NSET 

IR80 

IL«0 

CALL RANDDdX* IY*R0M) 

IXBIY 

X"Ta-XMlN ♦ ROM^XMN 
CALL PANI>D(IA«IH.«l)N) 

lAald 

YTaYMIN ♦ RDN»YMN 
X (N*l ) aX ( I » 

Y (M* 1 > ay d ) 
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I 


•j 


C CALCULATE T^E NUMHEP OF iMTEPSeCTION ALONG THE X AXIS 
no 3cC Kaj 

IFTrr.EQ,Y(K) .ANO.AT.EO.X (K)) GvO TO ni'' 

IF(YT.£Q.Y(K) .ANr).YT.eo.y(K*n ) go to lo 
tF(YT.FO.Y(K) ) GO TO 2v 
IF(r(K) ,GT.YT,AM1).YT.GT.Y(K*1) ) GO TO 4C 
IF (YT.Gf.Y (K) ,ANO.Y (K*l ) ,GT,YT) GO TO 40 
GO TO 3C'' 

l> IF(X(K) ,LT.XT*ANO.X(t<^l) ,6T,XT) GO TO OlC 
IF'X(K) ,GT.XT.ANO.x<K + n ,LT*XT) GO TO 010 
IF(XT-X(K)) lltSl'Wi? 

U IOaI»^NN(hf) 

GO TO OGO 
12 ILsIL>NN(K) 

GO TO 300 

?v IF(X(K)-XT) l2t310fU 

HI*. XXaX(K)* (YT-Y(K) )*(X(K^U-X.<Kn/(Y(KHl)-Y(Kn 
IF(XX-XT) 42«ai0*^l 

Hi 

GO TO 300 
42 ILaIL«>l 
300 CONTIMUF 

C CHECK WHETHER THE RANOOH POINT IS FALLING WITHIN THE BOUNDARY 
IF((IR-IR/2*?).E0.ft.OR.(lL-lL/2*2).EO.J) GO TO 302 
31 j NAaNA^l 
GO TO s:-) 

3C2 NORsNRP^l 

SOU continue 

AFa NA/<FL0AT(NA4.NRWn 
A«EAaAF*(XMAX-XMIN) * ( YMAX-YMIN) *FACT 
WRITE(NW.7) AFfAREA 

7 FORMAT (l.)x«2lHRELATIVE AREA RATIO a*FR*6t24H TOTAL AREA IN ACRES 
I a«Fl«>.6//> 

RETURN 

FNO 


SUPROUTINF RANOOnx*IY«YFL) 

C GENERATE THE RANDOM NUMHEP 

lYsIX* 40 '^R 
IF(Iy)5*6,6 
S IY 8 lY> 838 R 607 »l 
F. YFLsIY 

YFLsYFL/«38;l«)07.0 

RETURN 

EN(» 

I 



